table_function <- function(data_frames,
                          outcomes,
                          treatments,
                          trust_covariates,
                          interaction_effects = NULL,
                          interaction_effects_name = NULL,
                          additional_covariates = NULL,
                          fixed_effects_name,
                          weights_name = NULL,
                          table_name,
                          return_object = FALSE) {
  
  n_models <- length(data_frames)
  
  if (length(interaction_effects) == 0) {
      
      models <- lapply(1:n_models, function(i)
        
        if (identical(data_frames[[i]], diplomacy_long)) {
          lm_robust(as.formula(
            if (length(additional_covariates) == 0 & length(trust_covariates) != 0) 
            { paste0(outcomes[i], " ~ ", treatments[i]," + factor(", trust_covariates[i], ")", "*factor(country_exp)") }
            else if (length(additional_covariates) != 0 & length(trust_covariates) != 0) 
            { paste0(outcomes[i], " ~ ", treatments[i], " + factor(", trust_covariates[i], ")", "*factor(country_exp) +",
                     paste(additional_covariates, collapse = " + ")) }
            else if (length(additional_covariates) == 0 & length(trust_covariates) == 0) 
            { paste0(outcomes[i], " ~ ", treatments[i]) }
            else if (length(additional_covariates) != 0 & length(trust_covariates) == 0) 
            { paste0(outcomes[i], " ~ ", treatments[i], " + ",
                     paste(additional_covariates, collapse = " + ")) }
          ),
          fixed_effects = get(fixed_effects_name),
          clusters = mergeid,
          se_type = "stata",
          weights = if (length(weights_name) == 0) {NULL} else if (length(weights_name) != 0) {get(weights_name)},
          data = data_frames[[i]])
        } else {
          lm_robust(as.formula(
            if (length(additional_covariates) == 0 & length(trust_covariates) != 0)
            { paste0(outcomes[i], " ~ ", treatments[i], " + factor(", trust_covariates[i], ")") }
            else if (length(additional_covariates) != 0 & length(trust_covariates) != 0)  
            { paste0(outcomes[i], " ~ ", treatments[i], " + factor(", trust_covariates[i], ") +", paste(additional_covariates, collapse = " + ")) }
            else if (length(additional_covariates) == 0 & length(trust_covariates) == 0)
            { paste0(outcomes[i], " ~ ", treatments[i]) }
            else if (length(additional_covariates) != 0 & length(trust_covariates) == 0)  
            { paste0(outcomes[i], " ~ ", treatments[i], " +", paste(additional_covariates, collapse = " + ")) }
          ),
          fixed_effects = get(fixed_effects_name),
          se_type = "stata",
          weights = if (length(weights_name) == 0) {NULL} else if (length(weights_name) != 0) {get(weights_name)},
          data = data_frames[[i]]) 
        }
      )
      
  } else if (length(interaction_effects) == 1) {
    
    models <- lapply(1:n_models, function(i)
      
      if (identical(data_frames[[i]], diplomacy_long)) {
        lm_robust(as.formula(
          if (length(additional_covariates) == 0 & length(trust_covariates) != 0)
          { paste0(outcomes[i], " ~ ", treatments[i], "*", interaction_effects[[1]][i],
                   " + factor(", trust_covariates[i], ")", "*factor(country_exp)") }
          else if (length(additional_covariates) != 0 & length(trust_covariates) != 0)
          { paste0(outcomes[i], " ~ ", treatments[i], "*", interaction_effects[[1]][i],
                   " + factor(", trust_covariates[i], ")", "*factor(country_exp) + ",
                   paste(additional_covariates, collapse = " + ")) }
          else if (length(additional_covariates) == 0 & length(trust_covariates) == 0)
          { paste0(outcomes[i], " ~ ", treatments[i], "*", interaction_effects[[1]][i]) }
          else if (length(additional_covariates) != 0 & length(trust_covariates) == 0)
          { paste0(outcomes[i], " ~ ", treatments[i], "*", interaction_effects[[1]][i], " + ",
                   paste(additional_covariates, collapse = " + ")) }
        ),
        fixed_effects = get(fixed_effects_name),
        clusters = mergeid,
        se_type = "stata",
        weights = if (length(weights_name) == 0) {NULL} else if (length(weights_name) != 0) {get(weights_name)},
        data = data_frames[[i]])
      } else {
        lm_robust(as.formula(
          if (length(additional_covariates) == 0 & length(trust_covariates) != 0)
          { paste0(outcomes[i], " ~ ", treatments[i], "*", interaction_effects[[1]][i], " + factor(", trust_covariates[i], ")") }
          else if (length(additional_covariates) != 0 & length(trust_covariates) != 0)
          { paste0(outcomes[i], " ~ ", treatments[i], "*", interaction_effects[[1]][i], " + factor(", trust_covariates[i], ") + ", paste(additional_covariates, collapse = " + ")) }
          else if (length(additional_covariates) == 0 & length(trust_covariates) == 0)
          { paste0(outcomes[i], " ~ ", treatments[i], "*", interaction_effects[[1]][i]) }
          else if (length(additional_covariates) != 0 & length(trust_covariates) == 0)
          { paste0(outcomes[i], " ~ ", treatments[i], "*", interaction_effects[[1]][i], " + ",
                   paste(additional_covariates, collapse = " + "))}
        ),
        fixed_effects = get(fixed_effects_name),
        se_type = "stata",
        weights = if (length(weights_name) == 0) {NULL} else if (length(weights_name) != 0) {get(weights_name)},
        data = data_frames[[i]]) 
      }
    )
    
  } else if (length(interaction_effects) == 2) {
    
    models <- lapply(1:n_models, function(i)
      
          lm_robust(as.formula(
            if (length(additional_covariates) == 0 & length(trust_covariates) != 0)
            { paste0(outcomes[i], " ~ ", treatments[i], "*", interaction_effects[[1]][i], " + ", treatments[i], "*", interaction_effects[[2]][i], " + factor(", trust_covariates[i], ")") }
            else if (length(additional_covariates) != 0 & length(trust_covariates) != 0)
            { paste0(outcomes[i], " ~ ", treatments[i], "*", interaction_effects[[1]][i], " + ", treatments[i], "*", interaction_effects[[2]][i], " + factor(", trust_covariates[i], ") + ", paste(additional_covariates, collapse = " + ")) }
            else if (length(additional_covariates) == 0 & length(trust_covariates) == 0)
            { paste0(outcomes[i], " ~ ", treatments[i], "*", interaction_effects[[1]][i], " + ", treatments[i], "*", interaction_effects[[2]][i]) }
            else if (length(additional_covariates) != 0 & length(trust_covariates) == 0)
            { paste0(outcomes[i], " ~ ", treatments[i], "*", interaction_effects[[1]][i], " + ", treatments[i], "*", interaction_effects[[2]][i], " + ", paste(additional_covariates, collapse = " + ")) }
          ),
          fixed_effects = get(fixed_effects_name),
          se_type = "stata",
          weights = if (length(weights_name) == 0) {NULL} else if (length(weights_name) != 0) {get(weights_name)},
          data = data_frames[[i]])
    )
  }
    
    names(models) <- names(outcomes)
    
    if (length(weights_name) == 0) {  
      outcome_stats <- lapply(1:n_models, function(i)
        filter(model.frame(models[[i]]), get(treatments[i]) == 0) %>%
          dplyr::summarize(Mean = round(mean(get(outcomes[i]), na.rm = TRUE),2),
                    SD = round(sd(get(outcomes[i]), na.rm = TRUE),2))) 
    } else if (length(weights_name) != 0) {
      outcome_stats <- lapply(1:n_models, function(i)
        filter(model.frame(models[[i]]), get(treatments[i]) == 0) %>%
          dplyr::summarize(Mean = round(weighted.mean(get(outcomes[i]), w = get("(weights)"), na.rm = TRUE),2),
                    SD = round(sqrt(wtd.var(get(outcomes[i]), w = get("(weights)"), na.rm = TRUE)),2)))
    }
    
    outcome_stats <- lapply(1:n_models, function(i)
      outcome_stats[[i]] %>%
        mutate(Range = paste(sort(unique(model.frame(models[[i]])[,outcomes[i]])), collapse = ",")) %>%
        mutate(Range = paste0("\\{", Range, "\\}"))
    )
    
    names(outcome_stats) <- names(outcomes)
    
    treat_mean <- as.numeric(sapply(1:n_models, function(i)
      model.frame(models[[i]]) %>%
        dplyr::summarize(mean = round(mean(get(treatments[i]), na.rm = TRUE),2))
    ))
    
    if (length(interaction_effects) != 0) {
      
      interaction_stats <-
        lapply(1:length(interaction_effects), function (i)
          lapply(1:n_models, function(j)
            model.frame(models[[j]]) %>%
              dplyr::summarize(
                Mean = round(mean(get(
                  interaction_effects[[i]][j]
                ), na.rm = TRUE), 2),
                SD = round(sd(get(
                  interaction_effects[[i]][j]
                ), na.rm = TRUE), 2),
                Min = round(min(get(
                  interaction_effects[[i]][j]
                ), na.rm = TRUE), 2),
                Max = round(max(get(
                  interaction_effects[[i]][j]
                ), na.rm = TRUE), 2)
              )))
      
    }
    
    if (length(interaction_effects) == 0) {
    
      extra_rows <- rbind(
        sapply(outcome_stats, function(x) x[,'Range']),
        sapply(outcome_stats, function(x) x[,'Mean']),
        sapply(outcome_stats, function(x) x[,'SD']),
        treat_mean
      )
      
      extra_rows <- as_tibble(extra_rows) %>%
        mutate(term = c("Outcome range",
                        "Control outcome mean",
                        "Control outcome std. dev.",
                        paste0(unique(names(treatments)), " mean"))) %>%
        relocate(term)
      
      attr(extra_rows, 'position') <- 4:7
      
      treatment_names <- treatments[!duplicated(treatments)]
      coef_map_list <- setNames(names(treatment_names), treatment_names)
      
    } else if (length(interaction_effects) == 1) {
      
      extra_rows <- rbind(
        paste0("[", sapply(interaction_stats[[1]], function(x) x[,"Min"]),
               ",", sapply(interaction_stats[[1]], function(x) x[,"Max"]), "]"),
        sapply(interaction_stats[[1]], function(x) x[,'Mean']),
        sapply(interaction_stats[[1]], function(x) x[,'SD']),
        sapply(outcome_stats, function(x) x[,'Range']),
        sapply(outcome_stats, function(x) x[,'Mean']),
        sapply(outcome_stats, function(x) x[,'SD']),
        treat_mean
      )
      
      extra_rows <- as_tibble(extra_rows) %>%
        mutate(term = c(paste0(interaction_effects_name[[1]], " range"),
                        paste0(interaction_effects_name[[1]], " mean"),
                        paste0(interaction_effects_name[[1]], " std. dev."),
                        "Outcome range",
                        "Control outcome mean",
                        "Control outcome std. dev.",
                        paste0(unique(names(treatments)), " mean"))) %>%
        relocate(term)
      
      coef_map_list <- c(names(treatments[!duplicated(treatments)]),
                         rep(paste0(names(treatments[!duplicated(treatments)]), " $\\times$ ", interaction_effects_name[[1]]), n_models))
      
      names(coef_map_list) <- c(treatments[!duplicated(treatments)],
                                paste0(treatments[!duplicated(treatments)], ":", interaction_effects[[1]]))
      
      coef_map_list <- coef_map_list[!duplicated(names(coef_map_list))]
      
      attr(extra_rows, 'position') <- 6:13
      
    } else if (length(interaction_effects) == 2) {
      
      extra_rows <- rbind(
        paste0("[", sapply(interaction_stats[[1]], function(x) x[,"Min"]),
               ",", sapply(interaction_stats[[1]], function(x) x[,"Max"]), "]"),
        sapply(interaction_stats[[1]], function(x) x[,'Mean']),
        sapply(interaction_stats[[1]], function(x) x[,'SD']),
        paste0("[", sapply(interaction_stats[[2]], function(x) x[,"Min"]),
               ",", sapply(interaction_stats[[2]], function(x) x[,"Max"]), "]"),
        sapply(interaction_stats[[2]], function(x) x[,'Mean']),
        sapply(interaction_stats[[2]], function(x) x[,'SD']),
        sapply(outcome_stats, function(x) x[,'Range']),
        sapply(outcome_stats, function(x) x[,'Mean']),
        sapply(outcome_stats, function(x) x[,'SD']),
        treat_mean
      )
      
      extra_rows <- as_tibble(extra_rows) %>%
        mutate(term = c(paste0(interaction_effects_name[[1]], " range"),
                        paste0(interaction_effects_name[[1]], " mean"),
                        paste0(interaction_effects_name[[1]], " std. dev."),
                        paste0(interaction_effects_name[[2]], " range"),
                        paste0(interaction_effects_name[[2]], " mean"),
                        paste0(interaction_effects_name[[2]], " std. dev."),
                        "Outcome range",
                        "Control outcome mean",
                        "Control outcome std. dev.",
                        paste0(unique(names(treatments)), " mean"))) %>%
        relocate(term)
      
      coef_map_list <- c(names(treatments[!duplicated(treatments)]),
                         rep(paste0(names(treatments[!duplicated(treatments)]), " $\\times$ ", interaction_effects_name[[1]]), n_models),
                         rep(paste0(names(treatments[!duplicated(treatments)]), " $\\times$ ", interaction_effects_name[[2]]), n_models))
      
      names(coef_map_list) <- c(treatments[!duplicated(treatments)],
                                paste0(treatments[!duplicated(treatments)], ":", interaction_effects[[1]]),
                                paste0(treatments[!duplicated(treatments)], ":", interaction_effects[[2]]))
      
      coef_map_list <- coef_map_list[!duplicated(names(coef_map_list))]
      
      attr(extra_rows, 'position') <- 8:17
      
    }
    
    extra_rows <- rbind(extra_rows,
                        c("Number of blocks", sapply(models, function(i)
                          length(i$felevels[[1]]))))
    
    attr(extra_rows, 'position') <- c(attr(extra_rows, 'position'), max(attr(extra_rows, 'position'))+1)
    
    f <- function(x) format(round(x, 2), big.mark=",")
    gof_mapping <- list(
      list("raw" = "r.squared", "clean" = "R$^2$", "fmt" = f),
      list("raw" = "nobs", "clean" = "Observations", "fmt" = f)
    )
    
    modelsummary(models,
                 coef_map = coef_map_list,
                 add_rows = extra_rows,
                 gof_omit = "se_type|Adj.",
                 gof_map = gof_mapping,
                 stars = c("*" = 0.1, "**" = 0.05, "***" = 0.01),
                 escape = FALSE,
                 output = paste0("Tables/", table_name, ".tex"))
    
    if (return_object == TRUE) {
      
      table_object <- list(models = models,
                           extra_rows = extra_rows)
      
      return(table_object)
      
    } 
  
}